function doseresponse()
%
% compute dose-response based on existing data.
% use E_C for dose
% measure average of first 600 seconds for baseline
% measure average of 1500-2000 seconds for response (drug effect).
%
doselist = [0.1, 0.2, 0.4, 1.0]; % doses for PDAc D-R (or we could get them from E_C via unique...)
c = getuserdata; % get control list.
hc = findobj('tag', 'CTL');
sel = get(hc, 'value');
base_N2 = NaN * zeros(length(sel), length(doselist)); % baseline accumulation
base_P3 = NaN * zeros(length(sel), length(doselist)); % baseline accumulation
resp_N2 = NaN * zeros(length(sel), length(doselist)); % response accumulation
resp_P3 = NaN * zeros(length(sel), length(doselist)); % response accumulation
base_ppf_N2 = NaN * zeros(length(sel), length(doselist)); % response accumulation
resp_ppf_N2 = NaN * zeros(length(sel), length(doselist)); % response accumulation
base_ppf_P3 = NaN * zeros(length(sel), length(doselist)); % response accumulation
resp_ppf_P3 = NaN * zeros(length(sel), length(doselist)); % response accumulation
nobs = zeros(length(doselist), 1); % number of observations for each
%
k = zeros(length(doselist), 1);
for i = sel
   dose = str2num(c(i).E_C);
   j = find(doselist == dose);
   fp = c(i).FP; % get the FP structure
   t = fp.time;
   indx1 = find(t < 600);
   k(j) = k(j) + 1;
   base_N2 = mean(fp.maxN2(indx1));
   base_ppf_N2 = mean(fp.maxN22(indx1)./fp.maxN2(indx1));
   base_P3 = mean(fp.maxP3(indx1));
   base_ppf_P3 = mean(fp.maxP32(indx1)./fp.maxP3(indx1));
	indx2 = find(t > 1500 & t < 2000);   
   resp_N2 = mean(fp.maxN2(indx2));
   resp_ppf_N2 = mean(fp.maxN22(indx2)./fp.maxN2(indx2));
   resp_P3 = mean(fp.maxP3(indx2));
   resp_ppf_P3 = mean(fp.maxP32(indx2)./fp.maxP3(indx2));
   effect_N2{k(j), j} = resp_N2/base_N2;
   effect_P3{k(j), j} = resp_P3/base_P3;
   effect_ppf_N2{k(j), j} = resp_ppf_N2/base_ppf_N2;
   effect_ppf_P3{k(j), j} = resp_ppf_P3/base_ppf_P3;
   nobs(j) = nobs(j) + 1;
end;
for i = 1:length(doselist)
   m_N2(i) = mean([effect_N2{:,i}]);
   s_N2(i) = std([effect_N2{:,i}]);
   n_N2(i) = length([effect_N2{:,i}]);
   m_P3(i) = mean([effect_P3{:,i}]);
   s_P3(i) = std([effect_P3{:,i}]);
   n_P3(i) = length([effect_P3{:,i}]);
   m_ppf_N2(i) = mean([effect_ppf_N2{:,i}]);
   s_ppf_N2(i) = std([effect_ppf_N2{:,i}]);
   n_ppf_N2(i) = length([effect_ppf_N2{:,i}]);
   m_ppf_P3(i) = mean([effect_ppf_P3{:,i}]);
   s_ppf_P3(i) = std([effect_ppf_P3{:,i}]);
   n_ppf_P3(i) = length([effect_ppf_P3{:,i}]);
   
end;
s_N2 = s_N2 ./ sqrt(n_N2);
s_P3 = s_P3 ./ sqrt(n_P3);
s_ppf_N2 = s_ppf_N2 ./ sqrt(n_N2);
s_ppf_P3 = s_ppf_P3 ./ sqrt(n_P3);

h= findobj('tag', 'doseresponsecurve');
if(isempty(h))
   h = figure('tag', 'doseresponsecurve');
else
   figure(h);
   clf;
end;
subplot('Position', [0.05, 0.4, 0.4, 0.35]);
hold on;
plot(doselist, m_N2, 'marker', 'o', 'markersize', 6, 'markerfacecolor', 'k', 'color', 'k');
plot(doselist, m_P3, 'marker', 'v', 'markersize',6, 'markerfacecolor', 'white', 'color', 'k');
set(gca, 'Xlim', [0.095, 1.10]);
set(gca, 'Xscale', 'log');
ha = errorbar_log(doselist, m_P3, s_P3);
set(ha, 'color', 'k');
ha = errorbar_log(doselist, m_N2, s_N2);
set(ha, 'color', 'k');
set(gca, 'Ylim', [0.75, 2.25]);
yticks = [0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25];
set(gca, 'YTick', yticks);
for i = 1:length(yticks)
   ytl{i} = sprintf('%5.2f', yticks(i));
end;
set(gca, 'YTickLabel', ytl);
xticks = [0.1 0.2 0.4 1.0];
set(gca, 'XTick', xticks);
for i = 1:length(xticks)
   xtl{i} = sprintf('%5.1f', xticks(i));
end;
set(gca, 'XTickLabel', xtl);
set(get(gca, 'title'), 'string', 'Amplitude');

subplot('Position', [0.52, 0.4, 0.4, 0.35]);
hold on;
plot(doselist, m_ppf_N2, 'marker', 'o', 'markersize', 6, 'markerfacecolor', 'k', 'color', 'k');
plot(doselist, m_ppf_P3, 'marker', 'v', 'markersize',6, 'markerfacecolor', 'white', 'color', 'k');
set(gca, 'Xlim', [0.095, 1.10]);
set(gca, 'Xscale', 'log');
ha = errorbar_log(doselist, m_ppf_P3, s_ppf_P3);
set(ha, 'color', 'k');
ha = errorbar_log(doselist, m_ppf_N2, s_ppf_N2);
set(ha, 'color', 'k');
set(gca, 'Ylim', [0.25, 1.25]);
yticks = [0.25 0.5, 0.75, 1.00, 1.25];
set(gca, 'YTick', yticks);
for i = 1:length(yticks)
   ytl{i} = sprintf('%5.2f', yticks(i));
end;
set(gca, 'YTickLabel', ytl);
xticks = [0.1 0.2 0.4 1.0];
set(gca, 'XTick', xticks);
for i = 1:length(xticks)
   xtl{i} = sprintf('%5.1f', xticks(i));
end;
set(gca, 'XTickLabel', xtl);
set(get(gca, 'title'), 'string', 'PPF');

format short
[m_N2; s_N2; n_N2]'

[m_P3; s_P3; n_P3]'

x=[effect_N2];
x'

y=effect_P3;
y'



